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1 .  Introduction 


The  accurate  simulation  of  flows  with  features  spanning  many  length  scales  presents  a  major 
challenge  to  Computational  Fluid  Dynamics.  This  issue  becomes  especially  critical  in  3D 
flows  over  complex  configurations,  or  even  flows  over  simple  configurations  at  extreme  flow 
conditions,  where  multiple  interacting  3D  features  must  be  correctly  represented. 

The  last  decade  has  seen  the  emergence  of  two  basic  approaches  to  surmounting  this 
fundamental  difficulty.  High  resolution  upwind  schemes  attempt  to  represent  shocks  and 
contact  discontinuities  with  as  few  cells  as  possible  [1*6],  while  adaptive  grid  techniques 
automatically  adjust  the  local  mesh  dimension  in  an  attempt  to  resolve  the  scale  of  the  local 
flow  physics  [5-9].  Therefore,  the  combination  of  a  high  resolution  upwind  algorithm  with 
unstructured  mesh  adaptation  suggests  a  natural  and  extremely  flexible  platform  for 
simulating  flow  fields  with  features  spanning  a  variety  of  disparate  length  scales.  The  cunent 
work  adopts  this  combined  approach  and  uses  adaptively  refined  hexahedral  meshes  to 
provide  a  foundation  for  a  second  order  upwind  algorithm  based  upon  the  "Upwind  TVD" 
technique  of  Harten  and  Yee  [10], 

Mesh  adaptation  through  directional  cell  division  forms  a  flexible  and  efficient  procedure  for 
resolving  finer  length  scales  within  the  flow  field  [8,6].  The  technique  avoids  the  limited 
nature  of  grid  redistribution  algorithms,  and  minimizes  computational  expense  since  points 
are  only  added  in  regions  of  the  simulation  which  require  further  enhancement.  This 
flexibility  requires  a  solver  capable  of  dividing  an  arbitrary  collection  of  cells  and 
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necessitates  fully  unstructured  data  storage.  Despite  the  use  of  hexahedra!  meshes,  and  often, 
structured  looking  starting  meshes,  the  .solver  must  be  free  of  assumptions  about  mesh 
topology,  and  must  treat  any  mesh  as  simply  a  collection  of  cells  with  some  explicitly  defined 
connectivity. 

Although  interest  is  currently  growing  in  truly  multidimensional  {11]  and  multidirectional 
[12]  upwind  techniques,  most  successful  and  robust  upwind  techniques  still  rely  upon 
successive  application  of  essentially  one  dimensional  operators  through  space  operator 
splitting.  On  triangular  or  tetrahedral  meshes,  the  implementation  of  higher  order  upwind 
operators  (with  their  large  discretization  stencils)  is  not  straightforward,  and  .several 
approaches  have  been  suggested  in  the  recent  literature  [13,14].  The  use  of  hexahedral  cells, 
however,  provides  a  clear  mapping  of  this  discretization  stencil  onto  the  domain,  while 
avoiding  the  ambiguities  associated  with  higher  order  upwind  stencils  on  triangular  meshes 
[6].  Additionally,  such  meshes  easily  provide  the  high  quality  regular  mesh  near  wall 
boundaries  that  is  necessary  for  accurate  evaluation  of  the  viscous  fluxes  in  the  Navier-Stokes 
equations. 

Nearly  all  unstructured  techniques  require  an  initial  mesh  to  begin  the  numerical  simulation. 
In  the  current  method,  this  mesh  typically  begins  as  a  single  or  multiple  block  structured 
mesh.  Thus,  a  good  quality  background  mesh  must  still  be  generated  by  conventional 
techniques  and  the  flexibility  afforded  by  the  Delauny,  advancing  front,  or  other  techniques 
available  for  generating  tetrahedral  meshes  is  not  currently  available  for  hexahedral  meshes. 
This  drawback  is,  however,  no  more  severe  than  it  is  for  structured  techniques,  and  is 
largely  offset  by  both  the  natural  mapping  of  the  upwind  operator  and  accuracy  advantage 
offered  by  the  use  of  hexahedral  cells  in  vi.scous  region.s. 
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tion  of  Method 


2.  Descrio 


The  numerical  method  relies  upon  Harten  and  Yee’s  ‘‘Upwind  TVD”  algorithm  [10]  to 
integrate  the  inviscid  fluxes  while  the  viscous  terms  use  central  difference  modeling.  This 
spatial  discretization  is  advanced  in  time  through  a  modified  Runge-Kutta  procedure.  An 
option  for  central  differencing  of  the  inviscid  fluxes  also  exists.  This  option  is  essentially  a 
nodC'based  version  of  the  well  known  central  scheme  by  Jameson  with  blended  second  and 
fourth  order  artificial  dissipation  [15].  The  central  difference  option  was  included  primarily 
to  permit  direct  comparison  of  discrete  solutions  from  upwind  and  central  difference 
simulations  on  identical  adapted  meshes. 

2.a  Governing  Equations 


The  Navier- Stokes  equations  describe  the  unsteady  flow  of  a  viscous  gas  and  may  be  written 
in  integral  form  for  a  general  3D  region  V. 


V  dv 


F  *nds 


{1} 


F^Fi-Fy 


{2} 


Here,  U  is  the  state  vector  of  conserved  variables,  and  Fis  the  complete  tensor  of  flux  density 
which  contains  both  viscous  and  inviscid  components.  dV  is  the  closed  boundary  of  V  with 
the  outward  facing  unit  normal  vector  n  =  [njt,  riy,  This  set  of  equations  is  closed 
through  the  assumption  of  a  calorically  perfect  gas,  and  adopting  the  Stokes  Hypothesis. 


Assuming  an  isotropic,  Newtonian  fluid  yields  a  symmetric  shear  stress  tensor.  Constant 
Prandtl  number  modeling  is  assumed  throughout  the  domain,  and  Sutherland's  Law  descrihes 
the  variation  of  viscosity  with  temperature. 


2.b  Spatial  Integration  of  Convective  Fluxes 


The  spatial  discretization  of  the  governing  equations  may  be  ci>nveniently  developed  by 
focusing  on  the  symbolic  form  of  the  equations  given  by  expre.s,sion  1 1 }.  Considering  a 
control  volume  fixed  in  time,  and  applying  the  mean  value  theorem,  this  equation  may  he  re 
cast  as 


{3} 


The  inviscid  and  viscous  integration  implied  by  eq.{3}  proceed."!  on  a  dual  mesh  formed  by 
connecting  the  geometric  centers  of  the  physical  mesh  cells  which  surround  each  node. 
Figure  1  shows  construction  of  this  dual  auxiliary  mesh  in  two  and  three  dimensions.  The 
state  vector  of  conserved  quantities  is  stored  at  nodal  locations.  Expression  {3}  may  be 
further  specialized  for  application  to  a  specific  polyhedral  control  volume,  V,,  with  planar 
faces,  and  may  then  be  approximated  by: 


ai/y, 

dt 


1 

k  e  dVi 


^  I  Fvt 
*  *  € 


s, 


dt 


-V.  S  0. 

k  e  (IV, 


+  ^  X  rvt-St 

'  k  e  dV, 


{4} 


where  k  denotes  the  face  of  volume  V,  andS*  =  I5x,  Sy,  is  the  surface  vector  of  face  k. 
The  first  term  on  the  R.H.vS.  of  eq.{4}  approximates  the  inviscid  .surface  integral,  while  the 
second  term  balances  the  viscous  fluxes  through  the  facc.s  of  the  cell  around  /, 
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Auxiliary  Cell  i 


Physical  Cells 


V  ^ 

Physical  Cells 


Auxiliary  Cel!  i 


Figure  1.  Formation  of  dual  mesh  of  Auxiliary  Celts  in  two  and  three  dimensions. 


Evaluating  the  summatiiui  t\)r  the  inviscid  Hux  balance  mquircs  computation  of  the  numerical 
flux.  Qk,  through  each  face  k.  Using  the  nutation  found  in  Figure  2,  an  upwind  Uirmulation  o( 
this  flux  may  be  expressed  as 

Here  (  )l  and  {  )/?  denote  conditions  to  the  left  and  right  of  face  k  and  9?  is  the  right 
eigenvector  matrix  of  the  flux  Jacohians  in  transformed  space.  Eq.{5}  .separates  Qk  mio  a 
summation  of  a  symmetric  term  constructed  like  a  cenirai  difference  operator  and  a  term 
which  adapts  this  stencil  in  acetirdancc  with  local  wave  propagation. 

Tw-o  formulations  of  the  flux  function.  0f;,  were  suggested  by  Yee  which  lead  to  .second- 
order  TVD  schemes  11],  The  present  method  is  the  less  dissipative  of  the  two,  and  relics 
upon  the  “modified  flux”  approach  of  Hartcn.  With  /  =  l,...,.'i.  the  l'^  component  of  <J>i  may 
be  written  as 

kI  =  1  qA,')(x|  +  -  'M;.;  +  la' 

with  a*  = .  t7/J  j-jj 

leads  to  Harten  and  Yoe's  "Upwind  TVD"  scheme  [lOf  In  eq.{7}  lik  would  be  the 
difference  of  the  characteristic  variables  from  (  1/?  to  (  )i  if  9^’  were  evaluated  rn  the  right 
and  left  states  individually.  Instead,  the  nonlinear  nature  of  the  governing  equations  forces  us 
to  evaluate  the  mvcr.se  of  the  nghi  eigenvector  matrix  at  face  k  at  an  intermediate  slate  and  wc 
employ  the  approximate  Riemann  solver  developed  by  Roc  to  define  this  stale.  Here,  are 
the  eigenvalues  of  the  transformed  flux  Jacohians.  If  the  wave  speed,  A^.  were  to  vanish,  the 
asymmetric  term  on  the  R.H..S.  of  eq.I.Sl  would  also  vanish  and  thus  violate  the  entropy 
condition  (2|.  To  prevent  this.  lA  i  in  eq.{6}  is  ihrcsholdod  near  zero  by  introducing  the 
functions  ^iz)  and  }(- . 
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nzv 


and 


Ul 

id  ><5 

ld<^ 

2S 

U1 

!d>^ 

k|<5 

2S 

iih  S«  1 


{H} 


with  <>'«  1 


{9] 


The  limiter  function,  g^,  appearing  in  eqs.  {6}  and  {9}  is  taken  from  Ref.  1 16).  Notice  that 
setting  this  limiter  identically  to  zero  degenerates  the  entire  method  to  Roe’s  first  order 
scheme.  The  implementation  of  [17]  provides  a  basis  for  the  current  node-based  formulation 
(see  also  [6]). 


Implementation  of  Convective  Discretization 

Adopting  this  node-based  implementation  results  in  a  convenient  approach  for  the 
construction  of  an  unstructured  adaptive  algorithm,  while  still  permitting  accurate  and 
straightforward  treatment  of  boundaries.  In  developing  this  algorithm,  the  inviscid  flux 
balance  of  eq.{3}  is  evaluated  using  midpoint  quadrature,  and  thus  the  numerical  flux, 
must  be  evaluated  at  the  center  of  each  auxiliary  cell  face. 

Figure  3  shows  the  stencil  required  by  eqs. {3}  and  {6}  to  form  the  complete  numerical  flux 
through  the  west  face  of  V,.  These  expressions  emphasize  that  the  first  order  =  0  )  upwind 
approximation  of  the  inviscid  flux  requires  only  ncare.st  neighbor  information.  However,  as 
shown  in  the  exploded  view  of  this  stencil,  the  complete  flux  function  on  the  west  face  of  this 
cell,  -required  for  a  second-order  approximation  -  depends  on  information  from  beyond 
the  nearest  neighbor.  This  observation  reemphasizes  the  fundamental  difficulty  in  the  design 
of  higher  order  upwind  methods  on  unstructured  meshes. 

This  dilemma  is  alleviated  through  reexamination  of  the  formulation  pre.sented  by  eqs.  {6-9}. 
These  expressions  divide  the  formation  of  the  complete  stencil  into  a  sequence  of  operations 


8 


•  •  « 

8l  =  gR  r  g/?(  O^w » Ofj^ 


,  gi?) 


Figure  3.  Difference  stencil  for  second-order  upwind  representation  of  the  flux  function  <^at  the  west  face  of 

cell  i. 
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within  which  no  single  step  requires  communication  beyond  the  nearest  neighbor.  By  first 
calculating  O*  across  all  the  faces  of  all  cells  on  the  auxiliary  mesh,  and  then  storing  these 
values  in  each  cell,  the  limiters  Jk  -  and  therefore  the  flux  function  ^  -  may  be  evaluated  by 
simply  interrogating  only  the  auxiliary  cell  immediately  to  the  west  of  cell  V,  and  V,  itself. 
Breaking  up  the  second-order  discrete  operator  into  these  two  successive  steps  forms  the 
complete  update  to  the  state  vector  at  i  with  only  nearest  neighbor  inquiries. 

This  multiple  pass  evaluation  of  the  upwind  difference  stencil  is  analogous  to  the  splitting  of 
the  fourth  order  dissipation  operator  in  an  unstructured  central  difference  scheme.  In  such 
schemes,  the  fourth  difference  stencil  is  too  large  to  be  computed  with  only  nearest  neighbor 
information,  so  it  is  usually  evaluated  by  taking  the  second  difference  of  a  set  of  second 
differences  which  were  precomputed  in  a  previous  sweep  [9].  Similarly,  the  current 
procedure  forms  the  limiter  J  as  a  function  of  the  difference  of  characteristic  variables  that 
were  precalculated  in  a  prior  step.  The  additional  memory  requirements  implied  by  this 
strategy  are  not  severe  as  will  be  shown  in  a  subsequent  section. 

2.C  Spatial  Integration  of  Viscous  Fluxes 

The  inviscid  integration  scheme  may  be  extended  to  the  full  Navier-Stokes  equations  by  a 
separate  discretization  of  the  viscous  components  of  the  flux  density  tensor.  As  suggested  by 
the  fonn  of  eq.{4),  the  complete  update  to  any  node,  dj;,  is  a  summation  of  invi.scid  and 
viscous  contributions. 

4{7,  =  UZ7.)/  +  Ut7,)v  (!()} 

The  viscous  discretization  uses  central  difference  modeling  and  the  selection  of  control 
volumes  is  similar  to  that  found  in  Ref.  [8J.  The  choice  of  central  differencing  allows  the 
scheme  to  be  written  extremely  compactly,  and  eventually  it  involves  only  nodc-to-cell  and 
cell-to-node  communication  while  organizing  all  gather-scatter  operations  so  that  they  may 
be  easily  grouped  for  rapid  processing.  The.sc  communication  issues  remain  important 
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because  the  current  implementation  uses  a  cell  based  data  structure  (like  that  in  [8]  and  {?]) 
and  each  physical  cell  can  address  only  the  nodes  at  its  vertices  directly.  Thus  an 
implementation  based  upon  cell-to-node  operations  is  especially  convenient. 

Its  worth  noting  that  this  basic  discretization  is  also  amenable  to  edge-based  data  structures. 
In  fact.  Ref.  [18]  describes  the  formulation  of  Hessian  and  Laplacian  operators  using  only 
edge-based  formulas.  For  the  hexahedral  based  meshes  discussed  here,  however,  this  strategy 
holds  little  obvious  advantage,  and  so  the  cell-based  formulas  were  used. 


Ref.  [19]  contains  a  complete  development  of  the  viscous  discretization  in  two  and  three 
dimensions.  The  intent  of  this  section  is  to  provide  an  outline  of  that  discretization  so  that  the 
important  implementational  issues  may  be  clearly  understood.  The  complete  viscous 
discretization  may  be  illustrated  by  considering  the  formation  of  the  model  second  derivative 
term  Uxx  at  node  i.  Recalling  the  previous  assumption  that  the  control  volume  V,  has  planar 
faces  and  applying  the  Gauss  divergence  theorem  on  the  surface  dVi  of  auxiliary  cell  Vi  in 
Fig.l  gives 


(Uxx  )i  =  i 

/SV( 


(ux)nxClS 


=  -M(UxkSF-(Ux}wSj^  +  (UxkSi'-(Ux)sSi  +  (Ux)FSf-(Ux)BS^] 

Yi  [11] 


The  jc-components  of  both  the  unit  normal  vector,  rix,  and  the  surface  vector,  Sx,  are  taken 
oriented  in  positive  coordinate  directions.  With  the  auxiliary  cells  constructed  as  described 
above,  each  edge  incident  upon  node  i  will  always  pierce  the  face  of  an  auxiliary  cell.  Thus, 
eq.(  1 1 }  contains  one  contribution  from  each  edge  terminating  at  i.  This  observation  becomes 
increasingly  important  when  integrating  more  general  polyhedral  control  volumes. 

Eq.{  1 1 }  makes  use  of  first  derivatives,  Ux,  at  the  N,  S,  E,  W,  F,  B  faces.  A  second  surface 
integral,  similar  to  this  expression,  provides  each  of  these  quantities  at  the  required  midface 
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location.  This  is  accomplished  by  constructing  a  secondary  set  ol  control  volumes  which 
surround  each  of  the  edges  incident  upon  node  i.  Figure  4  illustrates  the  sevondary  veil 
constructed  around  the  edge  piercing  the  E  face  of  auxiliary  ceil  r  in  2  and  3D. 


urixdS 


=  -  ux,s{‘-  +  ujM  - 


Ve  is  the  volume  of  the  secondary  cell.  Eq.(  12}  requires  the  .surface  vectors  for  the  East 
secondai7  cell,  and  these  are  con.structed  by  averaging  the  surface  vectors  from  the  faces  of 
the  auxiliary  cells  surrounding  nodes  i  and  j  (the  nodes  which  define  the  end  points  of  the 
edge  which  pierces  the  E  face  of  Vj). 

=  +  =  +  Sf)  =  ((ff  +  9) 

=  ^(Sf  +  Sf )  ^  t  ^  Is  "')v,  =  J-lsf'  +  I , ,  I 

Similar  constructions  around  the  remaining  edges  which  meet  at  /  form  the  remaining  first 
derivatives  needed  to  evaluate  {11}.  The  choice  of  this  set  of  .secondary  control  volumes, 
and  the  use  of  the  Divergence  Theorem  en.sures  that  the  secondary  control  volumes  are 
closed,  and  that  the  viscous  discretization  remains  free  stream  preserving  on  any  arbitrary 
mesh.  Furthermore,  this  discretization  is  completely  conservative  and  the  formulation 
recognizes  and  damps  odd-even  decoupling  on  stretched  or  unstrctched  meshes  [  19]. 

Implementation  of  Viscous  Discretization 

Unstructured  implementation  of  the  discretization  outlined  by  eqs.j  1 1 )  and  {12}  e.s.scntially 
reduces  to  rearranging  these  expre.ssions  to  match  the  data  structure.  Within  the  current  cell- 
based  structure,  the  operations  comprising  these  two  surface  integrals  arc  broken  up  into 
contributions  from  each  physical  mesh  cell.  Operations  am  then  pt^rformed  within  each  me.sh 
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Figure  4.  Construction  of  Seconday  Cell  for  the  calculation  of  first  derivatives  on  the  E  face  of  the  auxiliary 

cell  which  surrounds  node  i.  (Left;  2D,  Right:  3D). 
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cell,  and  contributions  are  distributed  to  each  iit'  the  cell  vertices  (nodes)  This  process  !\irnis 
the  complete  update  at  the  mesh  nodes.  In  visuali/ing  this,  note  that  the  first  derivatives 
required  by  eq.{ll}  are  located  on  the  faces  of  the  auxiliary  cclis.  These  locations 
correspond  to  mid-edges  of  the  physical  mesh.  Thus,  the  algorithm  sweeps  through  the 
physical  mesh  on  a  cell-by-cel!  hasi.s,  forming  the  viscou.s  lluxes  along  the  edge  tif  the 
physical  cells,  and  finally  balancing  the.se  lluxes  on  the  auxiliary  cells  to  form  the  nodal 
update.  The  actual  implementation  is  di.scu.s.sed  fully  in  1 1*^]  and  is  ba.sed  upon  the  2D  work 
of  [20].  Implicit  in  this  approach  is  the  fact  that  under  the  assumptions  of  a  Newtonian  fduid, 
the  viscous  flux  calculation  becomes  a  linear  combination  of  first  derivative  quantities. 

2.d  Eigenvalue  Scaling  for  High  A.spect  Ratio  Cells 

Near  wall  boundaries,  high  a.spect  ratio  cells  arc  commonly  u.scd  to  efficiently  re.solvc  the 
boundary  layer.  Since  the  artificial  di.ssipation  is  .scaled  i,sotropically  with  the  spectral  radius 
of  the  flux  Jacobian  matrices,  the  damping  in  the  wall  nortr'>I  direction  may  become 
excessive.  Ref.  (21]  propo.sed  a  .'^D  extension  to  the  2D  variable  scaling  propo.sed  by 
Martinelli  [22]  which  adjusts  levels  of  the  blended  2nd  and  4th  order  smoothing  u.sed  in 
conjunction  with  the  central  difference  option  for  the  convective  llux  discretization.  When  the 
TVD  convective  dissipation  is  chosen,  this  noni.soiropic  scaling  affects  only  the  waves  which 
are  subject  to  the  entropy  cutoff.  In  the  current  work,  all  waves  are  cut  in  Euler  simulations, 
while  in  viscous  .simulations,  only  the  nonlinear  («  ±  a)  waves  are  .subject  to  this  cutoff.  In 
such  cases,  this  scaling  takes  a  form  similar  to  that  found  in  [2,^]. 
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3.  Adaptation 

The  adaptation  algorithm  increases  the  resolution  of  the  discrete  solution  by  locally  reducing 
the  mesh  spacing.  The  algorithm  begins  by  scanning  a  preliminary  (coarse)  solution  for 
regions  of  interest.  It  then  enhances  the  mesh  in  those  regions  through  the  addition  of  new 
computational  nodes.  These  new  points  are  currently  added  through  cell  division  [7],  and  the 
position  of  all  mesh  nodes  may  then  be  modified  by  a  redistribution  step  which  provides 
more  smoothly  varying  grid  metrics  without  altering  mesh  connectivity. 

3.a  Feature  Detection 

The  crisp  shocks  afforded  by  the  TVD  integration  scheme  may  be  detected  by  undivided 
differences  of  several  parameters.  However,  in  problems  involving  high  speed  flow  where 
several  parameters  may  jump  by  orders  of  magnitude  through  a  shock  (e.g.,  static  pressure), 
there  exists  a  danger  of  overlooking  less  prominent  features  with  a  crude  detection  algorithm. 
Furthermore,  recent  studies  of  various  feature  detection  strategies  demonstrate  that  failure  to 
adapt  to  smooth,  less  prominent  features  may  lead  to  an  adaptive  procedure  which  actually 
converges  to  the  wrong  solution  [23].  For  these  reasons,  a  new  procedure  has  been  developed 
which  separates  shock  recognition  from  the  detection  of  other  flow  features. 

Several  other  factors  contribute  to  this  decision.  The  current  algorithm  captures  shocks  with 
only  two  or  three  cells  [6,17]  because  the  difference  stencil  can  change  structure  to  reflect  the 
different  nature  of  the  governing  equations  on  either  side  of  such  an  interface.  With  a 
characteristic  thickness  of  only  a  few  mean  free  paths,  a  shock  is  one  of  the  very  few  flow 
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structures  which  will  always  be  much  smaller  than  the  local  mesh  dimension  (in  a  continuum 
situation).  Thus,  while  the  adaptation  will  eventually  embed  several  computational  cells 
across  most  other  flow  features,  this  will  not  occur  at  shocks.  Such  reasoning  gives  both 
physical  and  numerical  motivation  for  treating  shock  detection  in  a  .special  manner. 

Shock  Detection 

A  normalized,  undivided  second  difference  of  pressure  provides  a  parameter  which  is  quite 
sensitive  to  both  strong  and  weak  shocks.  Such  a  quantity  recognizes  curvature  in  the 
pressure  profile  across  a  three  point  difference  stencil.  Thus,  it  is  well  suited  to  recognize  the 
step-like  profiles  associated  with  the  discrete  shocks  expected  from  the  TVD  scheme  [1,17]. 
Additionally,  it  performed  well  in  the  study  of  refinement  parameters  provided  in  [7].  This 
second  difference  is  then  normalized  by  a  weighted  average  of  local  pressures  in  order  to 
reduce  its  magnitude  at  strong  shocks,  and  increase  its  size  at  weak  ones. 

Reference  [6]  contains  a  complete  formulation  of  this  shock  refinement  parameter,  Rs,  and 
any  physical  cell  with  Rs  exceeding  some  threshold  (usually  0.05)  is  tagged  for  division  as  a 
‘‘shock  cell."  Although  somewhat  empirical  in  basis,  this  procedure  very  reliably  locates 
shocks  in  trans-  super-  and  hypersonic  flow.  Moreover,  the  overall  detection  algorithm  is 
quite  insensitive  to  the  precise  threshold  value  chosen.  The  primary  drawback  of  this 
parameter  lies  in  the  fact  that  it  may  mistakenly  tag  a  strong  expansion  fan  as  a  “shock.” 
Nevertheless,  most  features  which  contain  strong  nonlinearity  usually  warrant  detection,  and 
this  is  not  a  serious  shortfall. 

Smooth  Feature  Detection 

With  the  cells  near  strong  nonlinearities  succc.ssfully  identified,  only  the  remaining  cells  are 
then  rescanned  for  smoother,  less  prominent  features.  With  these  extremea  removed,  the 
remaining  field  may  be  considered  “smooth,”  and  features  may  be  located  on  a  statistical 
basis  as  in  Ref.  [20].  Undivided  differences  of  total  velocity  (scaled  as  in  (23))  and  density 
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provide  a  basis  for  locating  inviscid  features*  A  second  sweep  using  an  undivided  second 
difference  of  velocity  magnitude  is  used  to  ivKate  regions  of  rapidly  varying  shear  stress  to 
identify  poorly  resolved  viscous  layers  in  Navier-Stokes  simulations.  All  physical  cells 
containing  values  of  the  refinement  parameter  above  a  threshold  set  20  percent  above  the 
mean  are  tagged  for  division. 

Directional  Division 

The  current  adaptive  procedure  identifies  not  only  a  feature’s  location,  but  also  its  orientation. 
Many  flow  features  are  predominantly  one-dimensional  and  are  frequently  almost  aligned 
with  the  coarse  base  grid.  With  this  information,  the  hexahedral  cells  may  be  divided 
directionally  to  avoid  excessive  resolution  in  directions  where  resolution  requirements  are 
more  relaxed.  In  three  dimensions,  this  means  that  a  cell  may  be  divided  in  only  one  or  two 
directions.  Such  refinements  create  only  one  or  three  additional  cells,  rather  than  the  seven 
which  are  created  if  each  cell  always  divides  in  all  three  directions. 

As  the  detector  scans  the  field,  it  stores  refinement  parameters  for  each  physical  cell  in  aH 
three  computational  directions.  Figure  5  contains  a  sketch  of  an  adaptation  map  (in  both  2 
and  3D)  within  which  all  the  cells  in  a  two-  or  three-dimensional  domain  will  lie.  Using  the 
notation  on  the  figure,  each  cell  will  have  some  [/?^,  /?t/.  coordinate  and  will  plot  on  this 
map.  The  threshold,  7,  mentioned  in  the  preceding  section  appears  as  a  quarter  sphere  on  this 
figure,  and  all  cells  outside  it  become  marked  for  refinement.  Furthermore,  the  location  of 
the  cells  on  this  map  indicate  the  directions  in  which  the  cells  require  refinement.  The  rays 
which  delineate  different  regions  emanate  from  the  origin  and  are  currently  set  at  an  angle  of 
15°. 

To  avoid  the  formation  of  excessive  interface  regions  between  cells  at  different  levels  of 
refinement,  directional  adaptation  was  used  only  at  the  finest  level  of  cell  division  [20,6}. 
This  final  division  typically  creates  from  50  percent  (2D)  to  75  percent  (3D)  of  the  final 
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nodes  in  the  mesh,  and  directional  division  lypicalfy  reduces  the  total  number  of  new  nodes 
created  by  about  50  precent  in  both  2  and  3D. 

4.b  Interfaces 

An  essential  complication  of  hexahedral  ceh  divi.sion  is  the  appearance  of  interfaces  between 
different  levels  of  physical  cells.  Figure  6  shows  a  region  of  divided  cells  bordering  on  an 
undivided  region  in  2  and  3D.  For  simplicity,  the  ceils  have  been  divided  in  all  directions. 
Shown  at  the  left  of  this  figure  are  the  physical  cells,  and  on  the  right  arc  the  auxiliary  cells 
formed  by  connecting  the  centers  of  the  physical  cells.  In  2D,  this  construction  results  in  the 
formation  of  triangular  auxiliary  cells  around  hanging  nodes  like  B.  In  3D  midface  nodes 
(like  B)  are  surrounded  by  a  pyramidal  auxiliary  cell  and  midedge  nodes,  like  C,  are 
surrounded  by  prismatic  cells.  All  of  these  auxiliary  cells  may  be  adopted  into  the  framework 
by  permitting  a  degeneracy  in  the  North  surface  vector  {  5^  s  0)  of  auxiliary  cells  B  and  C 
and  integrating  the  resulting  cells  without  special  treatment. 

Since  the  node-based  formulation  proceeds  by  interacting  the  auxiliary  cells,  without  special 
care,  all  cells  whose  difference  stencil  crosses  the  interface  will  suffer  a  first  order  stretching 
error  in  the  calculation  of  derivatives  across  this  interface.  Focusing  on  the  2D  example  in  the 
upper  half  of  the  figure,  this  means  that  although  cells  like  A  and  P  have  complete  difference 
stencils,  in  regions  where  flow  properties  vary  faster  than  linearly,  these  dilferences  will  be 
contaminated  by  mesh  stretching.  Triangular,  prismatic,  or  pyramidal  cells  also  suffer  a 
stretching  error,  but,  more  importantly,  their  difference  stencil  is  left  incomplete.  For 
example,  cell  B  has  no  northern  neighbor  which  gives  rise  to  a  degeneracy  in  the  northern 
extension  of  the  difference  stencil  of  B.  This  deficiency  is  overcome  by  linearly  interpolating 
for  the  missing®  (eq.{7})  from  theOyV  on  the  north  faces  of  A  and  C  (in  the  2D  case).  The 
treatment  of  6  as  a  degenerate  cell  whose  northern  surface  vector  is  zero  retains  storage  space 
for  the  interpolated  (.IJcn)b-  Such  a  quantity  is  required  by  the  limiter  on  the  south  face  of  cell 
B  when  second-order  accuracy  is  sought. 
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idcd  cells  in  two  and  three  dimensions 


These  ideas  form  a  simple  and  effective  strategy  for  treating  the  interfaces  which  appear  in 
adaptive  hexahedral  meshes.  First,  the  centers  of  all  physical  cells  are  used  to  define  the 
auxiliary  mesh.  Then,  if  a  degenerate  cell  i,s  encountered  while  sweeping  to  calculate  the 
difference  of  characteristic  variables,  a*,  the  mi.ssing  Uk  stem  from  a  combination  of  the  an 
on  similar  faces  of  its  neighboring  cells.  This  treatment  permits  a  complete  operator  to  be 
formed  on  the  face  opposite  the  degeneracy,  and  after  filling  in  the  missing  Ok,  the  fiux 
balance  for  the  prismatic  or  pyramidal  cells  proceeds  without  special  handling.  This  simple 
treatment  is  both  conservative  and  robust.  It  suffers  only  from  the  same  stretching  error 
experienced  by  any  node  whose  stencil  is  similarly  stretched.  The  viscous  update  to  the 
interface  nodes  proceeds  upon  the  same  set  of  control  volumes. 
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Numerical  Investigations  and  Results 


The  presentation  of  numerical  results  is  organized  to  first  establish  the  validity  of  the  solving 
procedure  and  then  to  examine  global  aspects  of  the  adaptive  methodology.  The 
investigations  focus  on  topologically  simple  example  problems  which  are  designed  to  first 
establish  the  accuracy  of  the  inviscid  and  viscous  discretization,  and  then  to  examine  the 
ability  of  the  adaptation  to  correctly  locate  and  resolve  structures  in  complex  3D  flows  with 
interacting  viscous  and  inviscid  features.  Unless  specifically  noted,  all  numerical  re.sults 
presented  use  the  Upwind  TVD  formulation  for  the  numerical  flux. 

4.  a  Fundamental  Issues 

This  section  begins  with  a  brief  examination  of  shock  and  boundary  layer  resolution, 
preservation  of  monotonicity,  and  other  basic  issues  of  the  numerical  modeling.  The  first 
numerical  example  considers  inviscid  flow  over  a  2D  circular  cylinder  exposed  to  a  Mach 
8.03  cross  flow.  The  left  half  of  Figure  7  contains  the  final  adapted  mesh  (1,500  nodes),  and 
line  plots  of  local  Mach  number  are  included  to  the  right.  These  Mach  number  distributions 
follow  a  ray  tracing  outward  from  the  centerline,  and  along  another  line  inclined  36®  from 
vertical.  As  is  evident  from  the  grid,  directional  division  was  employed  at  the  finest  level  of 
refinement. 

All  cases  presented  in  this  paper  were  computed  using  the  limiter  function,  taken  from 

Ref.  [16](Eq.  3.510 
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Figure  7.  Final  adapted  mesh  (1500  nodes)  and  Mach  No.  distributions  along  rays  A  and  B  for  inviscid  2D 

flow  over  circular  cylinder. 
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which  behaves  smoothly  and  is  not  overly  compressive. 


This  example  demonstrates  that  the  method  has  successfully  retained  the  favorable  shock 
capturing  properties  associated  with  the  Upwind  TVD  scheme  on  structured  meshes  (see  for 
example  Refs.  [10,  16,  17  and  4]).  Even  off  centerline  -  where  the  shock  is  not  grid  aligned  - 
the  discrete  shock  is  contained  within  two  cells.  Furthermore,  the  Mach  number  distributions 
in  each  plot  cross  four  grid  interfaces,  and  the  conservative  treatment  described  in  the 
preceding  section  successfully  maintains  a  smooth,  nonoscillatory  demeanor  through  these 
interfaces.  Measuring  to  the  sonic  point,  the  shock  stand-off  distance  in  this  example  agrees 
with  the  published  value  to  within  one  grid  point  [24]. 


A  second  inviscid  example  examines  the  full  3D  formulation  and  highlights  the  utility  of 
directional  adaptation  in  three  dimensions.  Figure  8  provides  an  overview  of  a  case  in  which 
a  double  wedge  comer  flow  is  simulated  at  Moo=2.98  with  wedge  angles  Sj=^=9.49°.  The 
final  mesh  contains  -70,000  nodes  which  is  14  times  fewer  than  the  roughly  980,000  which 
would  be  required  to  provide  the  same  shock  resolution  on  a  globally  refined  mesh. 
Directional  division  was  again  permitted  at  the  finest  mesh  level,  and  the  mesh  plane  depicted 
above  the  grid  (at  x=0.8)  shows  that  the  cells  have  adapted  themselves  to  match  the 
orientation  of  the  wedge  and  corner  shocks.  This  plane  appears  again  in  Figure  9,  where 
contours  of  constant  Mach  number  are  superimposed  directly  upon  both  experimental  data 
[25]  and  a  previous  inviscid  numerical  simulation  [26|.  Additionally,  the  undisturbed  wedge 
flow  away  from  the  corner  interaction  allows  direct  comparison  with  inviscid  gasdynamic 
theory.  The  agreement  with  experiment,  theory,  and  prior  calculations  results  is  very  good, 
and  the  shocks  span  no  more  than  two  cells.  Figure  10  .shows  a  convergence  history  for  this 
example  showing  a  net  reduction  in  the  global  residual  of  approximately  9  orders  of 
magnitude. 
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Figure  9.  Mach  contours  in  plane  at  jr  =  0.8  superimposed  upon  shock  positions  from  inviscid  gasdynamic 
theory  [24]  experimental  data  [25]  and  previous  inviscid  calculation  [26]. 
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Figure  10.  Convergence  history  of  RMS  sum  of  all  state  vector  residuals. 
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A  subsonic  flat  plate  boundary  layer  provides  an  initial  assessment  of  the  behavior  of  the 
viscous  discretization.  This  qua.si-2D  flow  was  computed  at  two  levels  of  resolution  using 
both  TVD  and  central  difference  operamrs  for  and  Rei=,5(XH).  This  simulatiiiO  used 

an  entropy  cutoff  in  the  TVD  scheme  of  0.1  and  a  fourth  difference  dissipation  coefficient  of 
•/l28  in  the  central  difference  calculation.  The  domain  extended  2L  upstream,  along,  and 
above  the  plate.  Figure  1 1  displays  u-velocity  profiles  of  both  .schemes  with  5  and  1.1  points 
in  the  boundary  layer.  With  only  five  points  in  the  boundary  layer,  the  Upwind  TVD  method 
has  already  essentially  reproduced  the  Bla.sius  profile.  This  result  contrasts  sharply  with  that 
from  the  central  difference  scheme  at  this  level  of  re.solution.  With  1 1  points  in  the  boundary 
layer,  the  two  integrations  produce  e.ssentially  identical  results.  At  this  higher  level  of 
resolution,  the  central  difference  result  produces  almost  no  evidence  of  the  “viscous 
overshoot”  evident  in  the  five  point  ca.se. 

Figure  12  shows  skin  friction  development  along  the  flat  plate.  This  plot  reports  results  for 
both  the  central  and  TVD  discretizations  at  two  levels  of  re.soJution  and  compares  these  with 
the  Blasius  relation  (^).664(/?er)''^‘.  With  13  points  in  the  boundary  layer  {Rei  of  5,{XX)), 
the  theoretical  skin  friction  is  predicted  wel!  by  about  Re^  =  100,  or  2  percent  of  the 
normalized  plate  length  from  the  leading  edge. 

4.b  Three-Dimensional  Viscous  Flow 

With  the  basic  properties  of  the  algorithm  established.  Figure  13  opens  the  presentation  of  a 
more  complex  example  aimed  at  facilitating  an  evaluation  of  the  feature  detection  algorithm. 
Success  or  failure  of  the  adaptive  strategy  rests  on  the  ability  of  this  routine  to  locate  and 
resolve  important  structures  and  smooth  regions  in  the  flow.  Thus,  it  is  important  to  evaluate 
the  feature  detection  algorithm  on  reali.stically  ctimplex,  viscous  flows  at  high  re.solution. 

The  test  case  beginning  in  Figure  13  considers  Mach  1.2  fiow  over  the  65°  cropped  delta 
wing  tested  in  Ref.  [2^]  at  20°  incidence  angle.  At  the.se  conditions,  this  flow  is  characteri/ed 
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Figure  1 1. 


Flat  plate  boundary  layer  calculation.  Mco  =  0.5,  Rei  =5,000.  Comparison  of  u-velocity  profile 
with  Blasius  solution  at  various  levels  of  mesh  resolution. 
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Figure  12.  Flat  plate  boundary  layer  calculation.  Moo  -  0.5,  Rei  =  5,000.  Comparison  of  skin  friction  with 
Blasius  solution  at  various  levels  of  mesh  rc.soluiion. 


by  a  strong  steady  primary  vortex  which  separates  at  the  sharp  leading  edge.  Underneath  this 
vortex,  a  cross-flow  shock  develops  which  induces  separation  of  the  boundary  layer,  and 
“locates”  the  secondary  separation.  The  tertiary  separation  lies  between  the  the  cross-flow 
shock  and  sharp  leading  edge.  On  the  windward  side,  a  bow  shctek  develops  to  turn  the  flow 
around  the  wing.  This  combination  of  interacting  viscous  and  inviscid  feature*  makes  this  a 
discriminating  problem  for  evaluating  the  capabilities  of  the  feature  detection  algorithm  for 
realistic,  complex  3D  flows. 

The  laminar,  half-span  calculation  used  a  Reynolds  number  of  48{)(X)0  based  on  the  root 
chord,  while  the  experiment  was  conducted  at  /?ec=2.4-5.3xl(K>.  The  wing  in  the  calculation 
was  not  truncated  at  the  trailing  edge,  but  instead,  was  extended  downstream  to  the  outflow 
boundary  (located  at  approx.  .r=l.lc). 

Figure  13  contains  an  overall  view  of  the  computation.  The  wing  is  depicted  with  the  mesh 
on  the  starboard  and  density  contours  on  the  port  side  of  the  delta  wing.  These  contours  show 
the  footprints  of  the  primary  vortex,  cross-flow  shock,  secondary  and  tertiary  vortices  on  the 
wing  surface.  The  final  adapted  mesh  used  999,000  nodes,  placing  about  260  points 
chordwise  and  150  points  spanwise  (at  the  trailing  edge)  on  the  upper  surface  of  the  wing. 
Figure  14  shows  further  details  of  the  flow  as  it  passes  through  the  plane  located  at  the 
trailing  edge  through  contours  of  total  pressure  Jos.s.  The  figure  clearly  shows  the  adaptation 
pattern  responding  to  capture  the  shock  triggered  separation  underneath  the  primary  vortex, 
as  well  as  resolution  of  feeding  sheet  and  the  entrained  tip  vortex.  This  indicates  that  Ute 
feature  detector  successfully  located  these  weaker  features,  while  still  responding  to  the 
strong  vortical  structures  and  bow  shock.  Figure  15  contains  symmetry  plane  Mach  contours, 
and  a  view  of  the  corresponding  adapted  grid  while  Figure  16  compares  this  calculation  with 
wind  tunnel  results  from  Ref.  [27].  The  computed  surface  streamlines  reproduce  the  overall 
surface  shear  pattern  from  the  experimental  oil  flow  visualization  including  formation  of  the 
secondary  and  tertiary  separations.  The  tertiary  vortex  in  the  calculation  does  not  appear  to 
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Figure  13.  Navier-Stokes  solution  of  supersonic  flow  over  delta  wing  of  Ref  [27]  at  M«,  =  1.2.  Rec  =  480,000, 
a  =  20°  with  999,000  nodes. 
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Figure  14.  Total  pressure  loss  contours  at  Le.  {x  =  1.0)  showing  shock  induced  separation. 


form  quite  as  early  or  as  fully  as  in  the  experinieni,  and  this  is  believed  to  be  a  Reynolds 
number  effect.  Figure  !6  also  contains  a  direct  comparison  of  the  Cp  distribution  at  jr=(),8c 
with  experimental  measurements.  Starting  at  the  symmetry  plane,  the  line  plot  agrees  with 
the  data  underneath  the  primary,  through  cross-llow  shock  and  secondary  separation  point, 
and  begins  to  differ  only  in  the  location  of  the  tertiary  separation.  Despite  the  Reynolds 
number  difference,  this  compari.son  is  quite  rea.sonahle.  Nevertheless,  it  should  be  re¬ 
emphasized  that  this  example  was  completed  to  focus  on  the  performance  of  the  detection 
algorithm  and  was  not  conducted  as  a  rigon)us  attempt  to  reproduce  wind  tunnel  data. 

4.C  Processing  Efficiency  and  Storage 

The  3D  Navier-Stokes  solver  has  been  ported  to  Cray  X-MP,  Y-MP,  and  Cray  2  platforms. 
The  main  iteration  loop  is  completely  vectorized,  including  the  numerical  flux  calculation 
and  the  computation  of  vi.scous  Ouxes.  No  explicit  edge  or  cell  coloring  schemes  were 
employed  to  achieve  this  result.  Instead,  a  coloring  is  implied  by  sequential  sweeps  over 
corresponding  faces  of  all  cells  (i.e.,  all  N  faces,  all  S  faces,  etc.).  Some  of  the  loops  within 
the  viscous  flux  calculation  required  “unwinding”  by  hand  to  vectorize  these  routines,  but  no 
other  special  coding  was  required. 

After  loading  a  restart  file  and  preprocessing  all  geometric  calculations  (surface  vectors,  cell 
volumes,  etc,),  timing  tests  were  run  on  5(K)  iterations  of  the  999,000  node  65°  cropped  delta 
wing  presented  in  the  preceding  .set  of  figures.  This  example  contained  50,(KK)  boundary 
nodes  and  the  boundary  conditions  were  not  fully  vectorized.  With  the  TVD  option  chosen, 
the  Navier-Stokes  simulation  proceeded  at  7()fi,sec/node/itcration  on  single  processor  of  the 
X-MP.  Selecting  central  differencing  for  the  inviscid  fluxes  results  in  Navier-Stokes 
simulations  at  30|i.scc/nodc/iteration.  Careful  timings  were  not  completed  on  either  the  Cray 
2  or  the  Y-MP. 
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The  soiver  currently  stores  153  words/node  for  lull  upwind  Navicr-Stokes  simulaiions.  30 
words/node  of  this  storage  are  required  for  the  precalculation  of  the  difference  of 
characteristic  variables  across  all  cel!  faces  (as  discussed  in  Section  2),  and  thus,  the  central 
difference  scheme  requires  approximately  125words/node.  On  a  per-node  basis,  the  scheme 
requires  3-4  times  more  storage  than  an  efficiently  written  structured  solver.  Adaptive 
methods  hope  to  offset  this  by  using  fewer  nodes  to  obtain  discrete  sedutions  of  comparable 
accuracy.  The  3D  examples  presented  earlier  in  this  .section  require  roughly  an  order  of 
magnitude  fewer  grid  points  than  an  equivalently  re.solved  structured  mesh  solution. 
Additionally,  the  use  of  fewer  nodes  not  only  reduces  storage  requirements,  but  also 
decreases  the  time  required  per  lime  step,  since  fewer  nodes  must  be  integrated.  With  the 
exception  of  local  time  stepping,  no  acceleration  techniques  were  applied  to  the  basic  solver. 
A  host  of  such  techniques  exist  for  accelerating  Runge-Kutta  schemes,  and  the  next  section 
discusses  likely  candidates  for  incorporation. 
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5.  Summary  and  Future  Work 


A  second-order  upwind  method  for  solution  of  the  3D  Navier-Stokes  equations  was 
formulated  on  adaptively  embedded  meshes.  This  implementation  permits  local  adaptation 
through  cell  division  in  response  to  emerging  flow  features.  The  feature  detection  algorithm 
separates  the  detection  of  shocks  from  that  of  other  flow  phenomena  and  can  divide  cells 
directionally  in  response  to  local  resolution  requirements. 

Application  of  the  method  to  several  inviscid  and  viscous  test  problems  pointed  out  the 
ability  of  the  method  to  provide  accurate  solution  to  examples  with  high  resolution 
requirements.  These  examples  demonstrated  that  the  unstructured  method  retains  the  crisp, 
nonoscillatory  shock  representation  associated  with  TVD  schemes  on  structured  meshes. 
Simulation  of  laminar,  viscous  flows  highlighted  the  ability  of  the  full  Navier-Stokes  solver 
to  accurately  model  viscous  flow  by  comparison  with  both  theoretical  and  experimental 
results  in  both  two  and  three  dimensions. 

With  the  preliminary  developmental  work  on  this  solver  complete,  further  efforts  will  focus 
on  improving  Uie  flexibility  and  utility  of  the  method.  The  Runge-Kutla  time  stepping  and 
embedded  mesh  topology  make  multigrid  acceleration  a  natural  choice  for  increasing  the 
convergence  efficiency  of  the  method.  Additionally,  such  work  .should  examine  pos.sibilities 
for  domain  decomposition  and  implicit  temporal  integration.  Finally,  future  investigations 
will  also  explore  the  practicality  of  applying  the  method  to  complex  configurations  using 
block-structured  initial  meshes. 
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Nomenclature 


^  small  parameter  preventing  zero  wave  speed.  (5«1 

e  total  internal  energy  per  unit  mass 

A  eigenvalue  of  flux  Jacobian  matrix 

M,  V,  w  Cartesian  velocity  componenus 
p  local  static  pressure 

Rs  shock  refinement  parameter 

p  density 

T  threshold  for  cell  division 

Vi  volume  of  cell  i 

F  complete  flux  density  tensor 

g  limiter  function 

n  unit  normal  vector  {rxx,  ny, 

5jk  surface  vector  of  face  k,  f5jt,  Sy,  5zj| 

Qk,  vector  of  numerical  llux  through  face  k 
U  state  vector  [p,  pu,  pv,  p\\\  e\ ' 

limited  flux  function  at  face  k 

91  right  eigenvector  matrix  of  Hux  Jacobian  in  transformed  space 
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